/**
 Dssv.cxx

 \file
 Implementation of class Dssv.

 \author Thomas Burton 
 \date 3/13/12
 \copyright 2012 BNL. All rights reserved.
*/

#include "tppmc/Dssv.h"

#include <cmath>
#include <fstream>
#include <stdexcept>

#include "tppmc/Fint.h"
#include "tppmc/Particle.h"

namespace tppmc {

   Dssv::Dssv() {
      try {
         InitialiseGrids();
      } // try
      catch(std::exception& e) {
         std::cerr << "Exception in Dssv::Dssv(): " << e.what() << std::endl;
      } // catch
   }

   Dssv::~Dssv() {
   }

   bool
   Dssv::InitialiseGrids() {
      double pointsQ2[Dssv::smNGridPointsQ2] = {
         0.8E0, 1.0E0, 1.25E0, 1.5E0, 2.E0, 2.5E0, 4.0E0, 6.4E0, 1.0E1, 1.5E1,
         2.5E1, 4.0E1, 6.4E1, 1.0E2, 1.8E2, 3.2E2, 5.8E2, 1.0E3, 1.8E3, 3.2E3,
         5.8E3, 1.0E4, 1.8E4, 3.2E4, 5.8E4, 1.0E5, 1.8E5, 3.2E5, 5.8E5, 1.0E6
      };
      mGridPointsQ2.assign(pointsQ2, pointsQ2 + GetNGridPointsQ2());
      
      double pointsX[Dssv::smNGridPointsX] = {
         1.0E-5, 4.0E-5, 6.7E-5, 1.0E-4, 1.4E-4, 2.0E-4, 3.0E-4, 4.5E-4, 6.7E-4,
         1.0E-3, 1.4E-3, 2.0E-3, 3.0E-3, 4.5E-3, 6.7E-3, 1.0E-2, 1.4E-2, 2.0E-2,
         3.0E-2, 4.5E-2, 0.06, 0.08, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25,
         0.275, 0.3, 0.325, 0.35, 0.375, 0.4,  0.45, 0.5, 0.55, 0.6, 0.65, 0.7,
         0.75,  0.8,  0.85, 0.9, 0.95, 1.0
      };
      mGridPointsX.assign(pointsX, pointsX + GetNGridPointsX());
      
      // Create a Grid with dimensions N_Q2 x (N_X - 1) for each flavour:
      Grid grid(GetNGridPointsQ2(), Array(GetNGridPointsX() - 1));
      mPartonGrids[LHAPDF::UP] = grid;
      mPartonGrids[LHAPDF::DOWN] = grid;
      mPartonGrids[LHAPDF::UBAR] = grid;
      mPartonGrids[LHAPDF::DBAR] = grid;
      mPartonGrids[LHAPDF::SBAR] = grid;
      mPartonGrids[LHAPDF::GLUON] = grid;
      // Read the grid values from a file
      std::ifstream ifs("tables/DSSV-GRID.NLO");
      
      for(int m = 0; m < GetNGridPointsX() - 1; ++m) {
         for(int n = 0; n < GetNGridPointsQ2(); ++n) {
            ifs
            >> mPartonGrids[LHAPDF::UP].at(n).at(m)
            >> mPartonGrids[LHAPDF::DOWN].at(n).at(m)
            >> mPartonGrids[LHAPDF::UBAR].at(n).at(m)
            >> mPartonGrids[LHAPDF::DBAR].at(n).at(m)
            >> mPartonGrids[LHAPDF::SBAR].at(n).at(m)
            >> mPartonGrids[LHAPDF::GLUON].at(n).at(m);
         } // for
      } // for
      
      LHAPDF::Flavour flavours[6] = {
         LHAPDF::UP, LHAPDF::DOWN, LHAPDF::UBAR,
         LHAPDF::DBAR, LHAPDF::SBAR, LHAPDF::GLUON
      };
      
      // Create the grids used in interpolation.
      // Note order of x and Q2 switched with respect to mPartonGrids
      const int nPoints = smNGridPointsX * smNGridPointsQ2;
      //      for(int i(1); i <= 6; ++i) {
      for(int i(0); i < 6; ++i) {
         // We treat the 2D x/Q2 grid as a 1D array for compatibility
         // with the CERN interpolation algorithm. The points are arranged as:
         //    (x_0,Q2_0),(x_1,Q2_0),...,(x_N-1,Q2_0)
         //    (x_0,Q2_1),(x_1,Q2_1),...,(x_N-1,Q2_1)
         //    (x_0,Q2_N-1),(x_1,Q2_N-1),...,(x_N-1,Q2_N-1)
         mInterpolationGrids.insert(std::make_pair(flavours[i], Array(nPoints)));
      } // for
      
      // Calculate the interpolation grid values
      for(int i = 0; i < GetNGridPointsQ2(); ++i) {
         for(int j = 0; j < GetNGridPointsX() - 1; ++j) {
            double XB0 = mGridPointsX.at(j);
            double XB1 = 1.0 - XB0;
            // We want to treat the interpolation grids as 2D arrays, but the
            // interpolation algorithm takes a 1D array as its argument so we
            // create a 1D array and determine the corresponding '2D' index.
            // The index is calculated as i * Nj + j to match the fortran
            // 2D array shaping used by the interpolation algorithm.
            int index = i * 47 + j;
            
            mInterpolationGrids[LHAPDF::UP].at(index)
            = mPartonGrids[LHAPDF::UP].at(i).at(j)
            / (::pow(XB1, 3) * XB0);
            
            mInterpolationGrids[LHAPDF::DOWN].at(index)
            = mPartonGrids[LHAPDF::DOWN].at(i).at(j)
            / (::pow(XB1, 4) * XB0);
            
            mInterpolationGrids[LHAPDF::UBAR].at(index)
            = mPartonGrids[LHAPDF::UBAR].at(i).at(j)
            / (::pow(XB1, 8) * ::pow(XB0, 0.5));
            
            mInterpolationGrids[LHAPDF::DBAR].at(index)
            = mPartonGrids[LHAPDF::DBAR].at(i).at(j)
            / (::pow(XB1, 8) * ::pow(XB0, 0.5));
            
            mInterpolationGrids[LHAPDF::SBAR].at(index)
            = mPartonGrids[LHAPDF::SBAR].at(i).at(j)
            / (::pow(XB1, 8) * ::pow(XB0, 0.5));
            
            mInterpolationGrids[LHAPDF::GLUON].at(index)
            = mPartonGrids[LHAPDF::GLUON].at(i).at(j)
            / (::pow(XB1, 5) * ::pow(XB0, 0.5));
         } // for
         
         // The last (N-1) x point at each Q2 point is set to zero:
         for(int k = 0; k < 6; ++k) {
            //         for(int k = 1; k <= 6; ++k) {
            LHAPDF::Flavour f = flavours[k];
            mInterpolationGrids[f].at(i * 47 + 46) = 0.0;//[i * 47 + 46] = 0.0;
         } // for
      } // for
      
      // Store the grid dimensions, to pass to the interpolator
      mGridDimensions.resize(2);
      mGridDimensions.at(0) = GetNGridPointsX();
      mGridDimensions.at(1) = GetNGridPointsQ2();
      
      // Store the grid points in a single array
      for(int i = 0; i < GetNGridPointsX(); ++i) {
         mGridPointsXQ2.push_back(::log(mGridPointsX.at(i)));
      } // for
      for(int i = 0; i < GetNGridPointsQ2(); ++i) {
         mGridPointsXQ2.push_back(::log(mGridPointsQ2.at(i)));
      } // for
      return true;
   }

   double Dssv::Interpolate(double x, double Q2, int flavour) const {
      std::vector<double> XT;
      XT.push_back(::log( x));
      XT.push_back(::log(Q2));
      Fint<double> fint;
      std::map<int, Array>::const_iterator array;
      array = mInterpolationGrids.find(flavour);
      return fint(2,
                  XT,
                  mGridDimensions,
                  mGridPointsXQ2,
                  array->second);
   }

   void Dssv::Evaluate(const Parton& p, double scale,
                       PartonContent& values) const {
      Parton parton(p); // Copy so we can change flavour
      parton.type = 2;
      values[LHAPDF::UP] = Evaluate(parton, scale, LHAPDF::UP);
      parton.type = 1;
      values[LHAPDF::DOWN] = Evaluate(parton, scale, LHAPDF::DOWN);
      parton.type = -2;
      values[LHAPDF::UBAR] = Evaluate(parton, scale, LHAPDF::UBAR);
      parton.type = -1;
      values[LHAPDF::DBAR] = Evaluate(parton, scale, LHAPDF::DBAR);
      parton.type = -3;
      values[LHAPDF::SBAR] = Evaluate(parton, scale, LHAPDF::SBAR);
      parton.type = 21;
      values[LHAPDF::GLUON] = Evaluate(parton, scale, LHAPDF::GLUON);
      values[LHAPDF::STRANGE] = values[LHAPDF::SBAR];
   }

   double Dssv::Evaluate(const Parton& parton, double Q2,
                         LHAPDF::Flavour flavour) const {
      // Perform bound checks on the x and Q2 point:
      try {
         if(not Supports(parton)) {
            return 0.;
         } // if
         if(parton.x < GetXMin() or parton.x > GetXMax()) {
            throw std::domain_error("x outside range valid for DSSV");
         } // if
         if(Q2 < GetQ2Min() or Q2 > GetQ2Max()) {
            throw std::domain_error("Q2 outside range valid for DSSV");
         } // if
         // We use the same distribution for strange and anti-strange:
         if(LHAPDF::STRANGE == flavour) {
            flavour = LHAPDF::SBAR;
         } // if
         // Run the interpolation to find xg1 at this (x,Q2)
         double xg1 = Interpolate(parton.x, Q2, flavour);
         // Apply flavour-dependent multiplication factors:
         double factor(0.);
         switch(flavour) {
            case LHAPDF::UP:
               factor = ::pow(1.0 - parton.x, 3) * parton.x;
               break;
            case LHAPDF::DOWN:
               factor = ::pow(1.0 - parton.x, 4) * parton.x;
               break;
            case LHAPDF::UBAR: // Deliberate fall-through - same for all qbar
            case LHAPDF::DBAR: // Deliberate fall-through - same for all qbar
            case LHAPDF::SBAR:
               factor = ::pow(1.0 - parton.x, 8) * ::pow(parton.x, 0.5);
               break;
            case LHAPDF::GLUON:
               factor = ::pow(1.0 - parton.x, 5) * ::pow(parton.x, 0.5);
               break;
            default:
               break;
         } // switch(flavour)
         return xg1 * factor / parton.x;
      } // try
      catch(std::exception& e) {
         std::cerr << "Exception: " << e.what() << std::endl;
         return 0.;
      } // catch
   }

   bool Dssv::Supports(const Parton& parton) const {
      bool isSupported(false);
      switch(parton.Lhapdf()) {
         case LHAPDF::DOWN :
         case LHAPDF::UP :
         case LHAPDF::STRANGE :
         case LHAPDF::DBAR :
         case LHAPDF::UBAR :
         case LHAPDF::SBAR :
         case LHAPDF::GLUON :
            isSupported = true;
            break;
         default:
            break;
      } // switch(flavour)
      return isSupported;
   }

   std::string Dssv::Name() const {
      return "DSSV helicity distributions";
   }
} // namespace tppmc
